#TEST#
#' Detect cloud, cloud shadow and snow using the Fmask method
#'
#' The Fmask (Function of mask) is an automated cloud, cloud shadow and snow masking detect ion method for Landsat TM/ETM+ images developed by Zhu and Woodcock (2012). A Matlab code for the Fmask is maintained by Zhu (http://code.google.com/p/fmask/). This is translation of the Matlab code to R. 
#'
#' @name FmaskforR-package
#' @aliases FmaskforR
#' @docType package
#' @title Automated cloud, cloud shadow and snow masking detection method for Landsat TM/ETM+ images.
#' @author Joseph Henry \email{joseph.henry@@sydney.edu.au}, Willem Vervoort \email{willem.vervoort@@sydney.edu.au}, Denis Duro and Mark Adams
#' @references Zhu, Z. and Woodcock, C.E. 2012. Object-based cloud and cloud shadow detection in Landsat imagery. Remote Sensing of Environment 118, 83-94. ' \url{http://code.google.com/p/fmask/}
#' @keywords package, Fmask, Landsat, Cloud, Snow, Shadow

### CONDITION OF USE ###
# Joseph reserves the right to use this code or its derivatives in the creation of an R package.

### ALPHA VERSION ###

FmaskForR <- function (INPUT_BANDS, INPUT_BANDS_DN, run.prob.over.water, L7.SUN.AZIMUTH.TAN.RADIANS, L7.SUN.AZIMUTH, L7.SUN.ELEVATION.TAN.RADIANS,OUTPUT) {
	#' Produce Landsat reflectance and temperature rasters.
	#' @param INPUT_BANDS file paths of reflectance and temperature data (bands 1 to 7) (character)
	#' @param INPUT_BANDS_DN file paths of digital number data (bands 1 to 3) (character)
	#' @param run.prob.over.water choice,  "yes" or  "no", to run Fmask over water. JH didn't use this because no large water bodies in catchment.
	#' @param OUTPUT file path of output file (character).
	#' @param L7.SUN.AZIMUTH.TAN.RADIANS sun azimuth in radians derived from metadata file that comes with tile (numeric)
	#' @param L7.SUN.AZIMUTH sun azimuth in degrees derived from metadata file that comes with tile (numeric)
	#' @param L7.SUN.ELEVATION.TAN.RADIANS  sun elevation derived from metadata file that comes with tile (numeric)

	################Constants set for testing 
		# INPUT_BANDS <- DATA_TR_OUT
		# INPUT_BANDS_DN <- DATA_PATH_RAW_EXCL
		# OUTPUT <- "C:/Users/Uni/Desktop/fmask.tif"
		# run.prob.over.water <- "no"
	################
		L7.FILLVALUE <- 255
		# Radiances
			Band_1_rad <- raster (INPUT_BANDS [1]) # Band 1
			gc()
			Band_1_rad_vec <- getValues(Band_1_rad) # Vectors
			gc()
			B1.OUT.PIX <- which (Band_1_rad_vec == Band_1_rad@data@min)
			rm (Band_1_rad)
			rm (Band_1_rad_vec)

			Band_2_rad <- raster (INPUT_BANDS [2]) # Band 2
			gc()
			Band_2_rad_vec <- getValues(Band_2_rad) # Vectors		
			gc()
			B2.OUT.PIX <- unique (c (B1.OUT.PIX, which (Band_2_rad_vec == Band_2_rad@data@min)))
			rm (Band_2_rad)
			rm (B1.OUT.PIX )
			rm (Band_2_rad_vec)
			gc()

			Band_3_rad <- raster (INPUT_BANDS [3]) # Band 3
			gc()
			Band_3_rad_vec <- getValues(Band_3_rad) # Vectors		
			gc()
			B3.OUT.PIX <- unique (c (B2.OUT.PIX, which (Band_3_rad_vec == Band_3_rad@data@min)))
			rm (Band_3_rad)
			rm (B2.OUT.PIX )
			rm (Band_3_rad_vec)
			gc()

			Band_4_rad <- raster (INPUT_BANDS [4]) # Band 4
			gc()
			Band_4_rad_vec <- getValues(Band_4_rad) # Vectors	
			gc()
			B4.OUT.PIX <- unique (c (B3.OUT.PIX, which (Band_4_rad_vec == Band_4_rad@data@min)))
			rm (Band_4_rad)
			rm (B3.OUT.PIX )
			rm (Band_4_rad_vec)
			gc()

			Band_5_rad <- raster (INPUT_BANDS [5]) # Band 5
			gc()
			Band_5_rad_vec <- getValues(Band_5_rad) # Vectors		
			gc()
			B5.OUT.PIX <- unique (c (B4.OUT.PIX, which (Band_5_rad_vec == Band_5_rad@data@min)))
			rm (Band_5_rad)
			rm (B4.OUT.PIX )
			rm (Band_5_rad_vec)
			gc()

			Band_6_rad <- raster (INPUT_BANDS [6]) # Band 6
			gc()
			Band_6_rad_vec <- getValues(Band_6_rad) # Vectors		
			gc()
			B6.OUT.PIX <- unique (c (B5.OUT.PIX, which (Band_6_rad_vec == Band_6_rad@data@min)))
			rm (Band_6_rad)
			rm (B5.OUT.PIX )
			rm (Band_6_rad_vec)
			gc()

			Band_7_rad <- raster (INPUT_BANDS [7]) # Band 7
			gc()
			Band_7_rad_vec <- getValues(Band_7_rad) # Vectors		
			gc()
			B7.OUT.PIX <- unique (c (B6.OUT.PIX, which (Band_7_rad_vec == Band_7_rad@data@min)))
			rm (Band_7_rad)
			rm (B6.OUT.PIX )
			rm (Band_7_rad_vec)
			gc()		
		
		# Digital numbers
			Band_1_DN <- raster (INPUT_BANDS_DN [1])[] # Band 1	
			B1Satu <- which (Band_1_DN == L7.FILLVALUE)  # Identify pixels that are saturated in the visible L7.PASS.CHOICE. Zhu has value of 255
			rm (Band_1_DN)
			gc()
			Band_2_DN <- raster (INPUT_BANDS_DN [2])[] # Band 2
			B2Satu <- which (Band_2_DN == L7.FILLVALUE)
			rm (Band_2_DN)
			gc()
			Band_3_DN <- raster (INPUT_BANDS_DN [3])[] # Band 3	
			B3Satu <- which (Band_3_DN == L7.FILLVALUE)	
			rm (Band_3_DN)
			gc()		
		## Clean
			B_Satu_unique_all <- sort (unique (c (B1Satu, B2Satu, B3Satu)))
			gc()
			B_Satu_table <- table( as.numeric (c (B1Satu, B2Satu, B3Satu)))
			dput (B1Satu, "B1Satu.RData"); dput (B2Satu, "B2Satu.RData"); dput (B3Satu, "B3Satu.RData");
			rm (B1Satu); rm(B2Satu); rm(B3Satu)
			gc()
			B_Satu_rows_all <- as.numeric (names (B_Satu_table[which (B_Satu_table == 3)]))
		 # Missing values
			id_missing <- B7.OUT.PIX 
			rm (B7.OUT.PIX)

		# Reflectance and temperature
			Band_1 <- raster (INPUT_BANDS [1]) # Band 1		
			gc()
			Band_1_vec_na <- Band_1[]
			rm (Band_1);  # Remove some of the rasters
			if (length (id_missing) > 0){ # Set NAs
				gc()
				Band_1_vec_na[id_missing] <- NA
			}
			dput(Band_1_vec_na, "Band_1_vec_na.RData")
			rm (Band_1_vec_na)

			Band_2 <- raster (INPUT_BANDS [2]) # Band 2		
			gc()
			Band_2_vec_na <- Band_2[]
			rm (Band_2);  # Remove some of the rasters
			if (length (id_missing) > 0){ # Set NAs
				gc()
				Band_2_vec_na[id_missing] <- NA
			}
			dput(Band_2_vec_na, "Band_2_vec_na.RData")
			rm (Band_2_vec_na)

			Band_3 <- raster (INPUT_BANDS [3]) # Band 3
			gc()
			Band_3_vec_na <- Band_3[]
			rm (Band_3);
			if (length (id_missing) > 0){ # Set NAs
				gc()
				Band_3_vec_na[id_missing] <- NA
			}
			dput(Band_3_vec_na, "Band_3_vec_na.RData")
			rm (Band_3_vec_na)
				  
			Band_4 <- raster (INPUT_BANDS [4]) # Band 4	
			gc()
			Band_4_vec_na <- Band_4[]
			rm (Band_4);	
			if (length (id_missing) > 0){ # Set NAs
				gc()
				Band_4_vec_na[id_missing] <- NA
			}	
			dput(Band_4_vec_na, "Band_4_vec_na.RData")
			rm (Band_4_vec_na)
				
			Band_5 <- raster (INPUT_BANDS [5]) # Band 5
			gc()
			Band_5_vec_na <- Band_5[]
			rm (Band_5);
			if (length (id_missing) > 0){ # Set NAs
				gc()
				Band_5_vec_na[id_missing] <- NA
			}	
			dput(Band_5_vec_na, "Band_5_vec_na.RData")
			rm (Band_5_vec_na)

			Band_6 <- raster (INPUT_BANDS [6]) # Band 6 in Celcius	
			gc()
			Band_6_vec_na <- Band_6[]
			rm (Band_6); 	
			if (length (id_missing) > 0){ # Set NAs
				gc()
				Band_6_vec_na[id_missing] <- NA
			}
			dput(Band_6_vec_na, "Band_6_vec_na.RData")
			rm (Band_6_vec_na)
				
			Band_7 <- raster (INPUT_BANDS [7]) # Band 7	
			gc()
			Band_7_vec_na <- Band_7[]
			rm (Band_7);	
			if (length (id_missing) > 0){ # Set NAs
				gc()
				Band_7_vec_na[id_missing] <- NA
			}
			dput(Band_7_vec_na, "Band_7_vec_na.RData")
			rm (Band_7_vec_na)
	print ("+++++++++++++++++++++++++++++++++++PART 1 +++++++++++++++++++++++++++++++++++")
	print(Sys.time())

	  ## Do		
		 # If all values are not NA
		  Band_7_vec_na <- dget( "Band_7_vec_na.RData")
		 if (length (id_missing) != length (Band_7_vec_na)){
			 rm (Band_7_vec_na)
			 Band_1_vec_na <- dget( "Band_1_vec_na.RData")
			 Band_2_vec_na <- dget( "Band_2_vec_na.RData")
			 Band_3_vec_na <- dget( "Band_3_vec_na.RData")
			 satu_Bv.data <- (Band_1_vec_na + Band_2_vec_na + Band_3_vec_na)[B_Satu_unique_all]
			 rm (Band_1_vec_na); rm (Band_2_vec_na); rm (Band_3_vec_na);
			 satu_Bv <- satu_Bv.data >= 1
			 rm ( satu_Bv.data) 
			 gc()
			 Satu_rows <- B_Satu_unique_all[which (satu_Bv == TRUE)]
			 rm (B_Satu_unique_all); rm (satu_Bv);
		# Potential cloud layer -- pass one
			Band_6_vec_na <- dget( "Band_6_vec_na.RData")
			BT <- Band_6_vec_na # Renamed (same as matlab name)
			rm (Band_6_vec_na);
				## Eq 1: Basic cloud test
					Band_3_vec_na <- dget( "Band_3_vec_na.RData");  Band_4_vec_na <- dget( "Band_4_vec_na.RData")
					NDVI <- (Band_4_vec_na - Band_3_vec_na) / (Band_4_vec_na + Band_3_vec_na) # Normalised Difference Vegetation Index
					dput(NDVI, "NDVI.RData")
					rm (Band_3_vec_na); rm (Band_4_vec_na); rm (NDVI);
					Band_2_vec_na <- dget( "Band_2_vec_na.RData");  Band_5_vec_na <- dget( "Band_5_vec_na.RData")
					NDSI <- (Band_2_vec_na - Band_5_vec_na) / (Band_2_vec_na + Band_5_vec_na) # Normalised Difference Snow Index	
					rm (Band_2_vec_na); rm (Band_5_vec_na);
					NDVI <- dget("NDVI.RData")		
					Band_7_vec_na <- dget( "Band_7_vec_na.RData")			
					Basic_Test <- ifelse (Band_7_vec_na > 300 & BT < 2700 & NDSI < 0.8 & NDVI < 0.8, TRUE, FALSE) # Run test
					dput (Basic_Test,"Basic_Test.RData");  rm (Basic_Test); rm (Band_7_vec_na); rm (BT)
					dput (NDSI,"NDSI.RData");  rm (NDSI);
					gc()
				## Eq 2: Whiteness test
					Band_1_vec_na <- dget( "Band_1_vec_na.RData"); Band_2_vec_na <- dget( "Band_2_vec_na.RData"); Band_3_vec_na <- dget( "Band_3_vec_na.RData");
					MeanVis <- (Band_1_vec_na + Band_2_vec_na + Band_3_vec_na) / 3
					B1.MeanVis.diff <- abs (Band_1_vec_na - MeanVis)
					rm (Band_1_vec_na)
					B2.MeanVis.diff <- abs (Band_2_vec_na - MeanVis)
					rm (Band_2_vec_na)
					B3.MeanVis.diff <- abs (Band_3_vec_na - MeanVis)
					rm (Band_3_vec_na)
					Bsum.MeanVis.diff <- (B1.MeanVis.diff + B2.MeanVis.diff + B3.MeanVis.diff) / MeanVis
					rm (MeanVis); rm (B1.MeanVis.diff); rm (B2.MeanVis.diff); rm (B3.MeanVis.diff);
					whiteness <- Bsum.MeanVis.diff
					rm (Bsum.MeanVis.diff)
					gc()
				## Adjust whiteness values
					# Set saturated whiteness to zero
					if (length (Satu_rows) > 0){
						whiteness[Satu_rows] <- 0
					}
					whiteness[whiteness > 1] <- 1 #C- whiteness goes over 1, but I don't know why. Therefore reset to 1. -C#	
					Whiteness_Test <- whiteness < 0.7
					rm(whiteness);
					dput(Whiteness_Test,"Whiteness_Test.RData")
					gc()
				## Eq 3: Haze test
					Band_1_vec_na <- dget( "Band_1_vec_na.RData"); Band_3_vec_na <- dget( "Band_3_vec_na.RData");
					HOT_Test.data <- Band_1_vec_na - 0.5 * Band_3_vec_na - 800 # Test values
					HOT_Test <- HOT_Test.data > 0 # Run test
					rm (Band_1_vec_na); rm (Band_3_vec_na); rm (HOT_Test.data)
					dput(HOT_Test,"HOT_Test.RData")
					rm (HOT_Test)
					gc()
				## Eq 4: Band 4 and 5 ratio
					Band_4_vec_na <- dget( "Band_4_vec_na.RData"); Band_5_vec_na <- dget( "Band_5_vec_na.RData");
					B45_Test <- (Band_4_vec_na / Band_5_vec_na) > 0.75 # Test
					rm (Band_4_vec_na); rm (Band_5_vec_na)
					dput (B45_Test, "B45_Test.RData")
					rm (B45_Test)
					gc()
				## Eq 5: Water test
					NDVI <- dget("NDVI.Rdata"); Band_4_vec_na <- dget( "Band_4_vec_na.RData");
					Water_Test <- ifelse ((NDVI < 0.01 & Band_4_vec_na < 1100)|(NDVI < 0.1 & NDVI > 0 & Band_4_vec_na < 500), TRUE, FALSE)
					rm (NDVI); rm (Band_4_vec_na);
					dput (Water_Test, "Water_Test.RData")
					rm (Water_Test)
					gc()
				## Eq 6: Potential cloud ?
					Basic_Test <- dget("Basic_Test.RData"); Whiteness_Test <- dget("Whiteness_Test.RData"); HOT_Test <- dget("HOT_Test.RData"); B45_Test <- dget("B45_Test.RData");
					PCP <- ifelse (Basic_Test == TRUE & Whiteness_Test == TRUE & HOT_Test == TRUE & B45_Test ==TRUE, TRUE, FALSE)
					rm (Basic_Test); rm (Whiteness_Test); rm (HOT_Test); rm (B45_Test);		
					gc()
				
				####S Choice switch
					choice_switch_1 <- length (which (PCP == TRUE)) / length (PCP)
					dput(PCP,"PCP.RData")
					rm (PCP)
					gc()
				####E
				
				if (choice_switch_1 >= 0.99) { # Too many clouds, go straight to end. Set entire tile to cloud.
					Band_1_vec_na <- dget( "Band_1_vec_na.RData"); 
					Band_1 <- raster (INPUT_BANDS [1]);
					Fmask <- rep (1, length (Band_1_vec_na)) # Set vector
					Fmask_ras <- setValues (Band_1, Fmask) # Turn vector into raster
					writeRaster(Fmask_ras,OUTPUT,overwrite = TRUE) # Save raster
					rm (Band_1_vec_na); rm (Band_1); rm (Fmask);rm (Fmask_ras);
				}	
	print ("+++++++++++++++++++++++++++++++++++PART 2 +++++++++++++++++++++++++++++++++++")
	print(Sys.time())

				if (choice_switch_1 < 0.99){ # Not fully cloud. Potential cloud layer -- pass two 
					# Temperature probability for water
					if (run.prob.over.water == "yes"){ # Are there any bodies of water to run this over? 
						## Eq 7: Temperature test over water
							Water_Test <- dget("Water_Test.RData")
							Band_7_vec_na <- dget ("Band_7_vec_na.RData")
							Clearsky_Water <- (Water_Test == TRUE & Band_7_vec_na <= 300) #C- Inconsistency between matlab and paper? -C#
							rm (Water_Test); rm (Band_7_vec_na);
						## Eq 8
							Band_6_vec_na <- dget( "Band_6_vec_na.RData")
							BT <- Band_6_vec_na # Renamed (same as matlab name)
							rm (Band_6_vec_na);
							T_water <- quantile (as.numeric (BT[Clearsky_Water == TRUE]), 0.825, na.rm = TRUE)
							dput (T_water,"T_water.RData")
							gc()
						## Eq 9
							wTemperature_Prob <- (T_water - BT) / 400
							rm (BT)
							dput (wTemperature_Prob, "wTemperature_Prob.RData")
							rm (wTemperature_Prob)
							rm (T_water)
							gc ()
						## Eq 10: Brightness test over water
							t_bright <- 1100
							Band_5_vec_na <- dget ("Band_5_vec_na.RData")
							Brightness_Prob <- Band_5_vec_na / t_bright
							rm (Band_5_vec_na)
							Brightness_Prob[Brightness_Prob > 1] <- 1
						## Eq 11: Final probability mask over water
							T_water <- dget ("T_water.RData")
							if (is.na (T_water) == FALSE){
								rm (T_water)
								wTemperature_Prob <- dget ("wTemperature_Prob.RData")
								wCloud_Prob <- wTemperature_Prob * Brightness_Prob * 100
								rm (Brightness_Prob); rm (wTemperature_Prob)
								gc ()
								} else { #C- If no clear-sky water was identified -C#
								rm (T_water)
								wCloud_Prob <- 1 / (Brightness_Prob) 
								rm (Brightness_Prob)
								dput (wCloud_Prob,"wCloud_Prob.RData")
								gc()
							}
						}
				# Temperature probability for land
					## Eq 12
						T_buffer <- 4 * 100
						PCP <- dget ("PCP.RData"); Water_Test <- dget ("Water_Test.RData")
						Clearsky_Land <- ifelse (PCP == FALSE & Water_Test == FALSE, TRUE, FALSE)
						rm (PCP); rm (Water_Test);
						dput(Clearsky_Land,"Clearsky_Land.RData")
						gc()
					## Eq 13	
						Band_6_vec_na <- dget( "Band_6_vec_na.RData")
						BT <- Band_6_vec_na # Renamed (same as matlab name)
						rm (Band_6_vec_na);
						T_low <- quantile (as.numeric (na.omit (BT[Clearsky_Land == TRUE])), 0.175) # Low percentile
						T_high <- quantile (as.numeric (na.omit (BT[Clearsky_Land == TRUE])), 0.825) # High percentile
						rm (Clearsky_Land); rm (BT)
					# Other stuff
						t_tempL <- T_low - T_buffer
						t_tempH <- T_high + T_buffer
						rm (T_high); 
						dput(T_low, "T_low.RData"); rm (T_low);
						Temp_l <- t_tempH - t_tempL
						dput(t_tempL,"t_tempL.RData"); rm (t_tempL)
						
					## Eq 14
						Band_6_vec_na <- dget( "Band_6_vec_na.RData")
						BT <- Band_6_vec_na # Renamed (same as matlab name)
						rm (Band_6_vec_na);
						lTemperature_Prob <- (t_tempH - BT) / Temp_l
						rm (BT)
						dput(t_tempH,"t_tempH.RData"); rm (t_tempH)
						dput(Temp_l,"Temp_l.RData"); rm (Temp_l)
						lTemperature_Prob[lTemperature_Prob < 0] <- 0 #C- Probability can be > 1 -C#
					####S Choice switch
						Clearsky_Land <- dget("Clearsky_Land.RData")
						choice_switch_2 <- length (which (Clearsky_Land == TRUE)) / length (Clearsky_Land)
						rm (Clearsky_Land)
					####E
				}

					#C- Fmask will use the clear clear-sky pixels from both land and water for calculating temperature probability -C#
					if (choice_switch_1 < 0.99){
						if (choice_switch_2 < 0.001){ 	
						# NDSI
							NDSI <- dget ("NDSI.RData")
							modified_NDSI <- NDSI
							B3Satu <- dget ("B3Satu.RData")
							modified_NDSI_tf_1 <- NDSI[B3Satu]<0
							rm (NDSI)
							modified_NDSI_tf_2 <- B3Satu[which(modified_NDSI_tf_1==TRUE)]
							rm (B3Satu); rm (modified_NDSI_tf_1)
							modified_NDSI[modified_NDSI_tf_2] <- 0								
						## Eq 20: potential snow layer
							Band_4_vec_na <- dget ("Band_4_vec_na.RData"); Band_2_vec_na <- dget ("Band_2_vec_na.RData"); 
							Band_6_vec_na <- dget( "Band_6_vec_na.RData")
	print ("+++++++++++++++++++++++++++++++++++PART 3 +++++++++++++++++++++++++++++++++++")
	print(Sys.time())

							BT <- Band_6_vec_na # Renamed (same as matlab name)
							rm (Band_6_vec_na);
							PSL <- ifelse(modified_NDSI>0.15&BT<380&Band_4_vec_na>1100&Band_2_vec_na>1000,TRUE,FALSE)
						# Vector with values 4) cloud pixels, 2) shadow pixels, 3) snow pixels
							Band_1 <- raster (INPUT_BANDS [1])
							Fmask <- rep(0,length(Band_1))
							rm (modified_NDSI); rm (Band_4_vec_na); rm (Band_2_vec_na); rm (BT);
							Fmask[PSL == TRUE] <- 4 	# Snow
							Fmask[id_missing] <- NA # NA post 2003
							Fmask_ras <- setValues(Band_1,Fmask) # Create raster
							writeRaster(Fmask_ras,OUTPUT,overwrite = TRUE) # Save raster
							rm (Band_1); rm (Fmask_ras);
						}
					}

				#C- Fmask will use the clear clear-sky pixels land only for calculating temperature probability -C#
					if (choice_switch_1 < 0.99){				
						if (choice_switch_2 >= 0.001){			
							## Variability probability
							 # Modify vales for saturated pixels
							  # NDVI
								NDVI <- dget ("NDVI.RData")
								modified_NDVI <- NDVI	
								B2Satu <- dget ("B2Satu.RData")
								modified_NDVI_tf_1 <- NDVI[B2Satu] > 0
								modified_NDVI_tf_2 <- B2Satu[which (modified_NDVI_tf_1 == TRUE)]
								rm(modified_NDVI_tf_1); rm (B2Satu) ; rm (NDVI)
								modified_NDVI[modified_NDVI_tf_2] <- 0
								dput(modified_NDVI,"modified_NDVI.RData")
								rm (modified_NDVI)
								# NDSI
									NDSI <- dget ("NDSI.RData")	
									modified_NDSI <- NDSI
									B3Satu <- dget ("B3Satu.RData")
									modified_NDSI_tf_1 <- NDSI[B3Satu] < 0
									rm (NDSI); 
									modified_NDSI_tf_2 <- B3Satu[which (modified_NDSI_tf_1==TRUE)]
									rm (B3Satu); rm (modified_NDSI_tf_1);
									modified_NDSI[modified_NDSI_tf_2] <- 0
									rm (modified_NDSI_tf_2);
									dput (modified_NDSI,"modified_NDSI.RData")
									rm (modified_NDSI);
								## Eq 15
	#START###################################################
		# COMMENT
		# JH (25/4/12): Change from the Matlab or paper.					
								#C- I removed whiteness because it can equal 1. Hence, Variability prob would be 0. Therefore, subsequent probabilities would be 0 -C#
									modified_NDVI <- dget("modified_NDVI.RData"); modified_NDSI <- dget ("modified_NDSI.RData")
									Variability_Prob <- 1 #-max(abs(modified_NDVI),abs(modified_NDSI),na.rm=T)#,whiteness,na.rm=T)
									rm (modified_NDSI); rm (modified_NDVI)
									#Variability_Prob <- 1-max(abs(modified_NDVI))#,abs(modified_NDSI))
	#END#####################################################							
							## Cloud probability for land
								## Eq 16
									lCloud_Prob <- (lTemperature_Prob * Variability_Prob) * 100
									rm (Variability_Prob)
								 ## Eq 17
									cldprob <- 22.5 # Default cloud probability threshold for cloud detection
									Clearsky_Land <- dget ("Clearsky_Land.RData")
									Land_threshold <- as.numeric (quantile (lCloud_Prob[Clearsky_Land == TRUE], 0.825, na.rm = TRUE) + cldprob) 
									rm (Clearsky_Land);
									# Treshold for defining cloud cover
						  ## Potential cloud layer
							 ## Eq 18
								PCP <- dget ("PCP.RData");  Water_Test <- dget ("Water_Test.RData");
								if (run.prob.over.water=="yes"){							
									clr_max <- 50 # Treshold
									wCloud_Prob <- dget ("wCloud_Prob.RData")
									PCL_t_1 <- PCP == TRUE & Water_Test == TRUE & wCloud_Prob > clr_max # Cloud over water
									rm (wCloud_Prob)
									PCL_t_1 <- dput (PCL_t_1, "PCL_t_1.RData")
									rm (PCL_t_1)
								}
								PCL_t_2 <- PCP == TRUE & Water_Test == FALSE & lCloud_Prob > Land_threshold # Cloud over land
								rm (Land_threshold); rm (PCP);
								PCL_t_2.length <- length (PCL_t_2)
								PCL_t_2 <- dput (PCL_t_2, "PCL_t_2.RData")
								rm (PCL_t_2)
								PCL_t_3 <- lCloud_Prob > 99 & Water_Test == FALSE # High probability of cloud over land
								rm (Water_Test); rm (lCloud_Prob);
								PCL_t_3 <- dput (PCL_t_3, "PCL_t_3.RData")
								rm (PCL_t_3)
								BT <- dget( "Band_6_vec_na.RData") # Band6 and BT are the same thing, just change of nomenclature
								T_low <- dget ("T_low.RData")
								PCL_t_4 <- BT < (T_low - 3500) # Really cold cloud
								rm (BT); rm (T_low);
								PCL_t_4 <- dput (PCL_t_4, "PCL_t_4.RData")
								rm (PCL_t_4)
							## Classify each pixel
								PCL_raw_class <- rep (NA, PCL_t_2.length) # Blank vector
								rm (PCL_t_2.length)
								# Add values to vector
									if (run.prob.over.water=="yes"){
										PCL_t_1 <- dget ("PCL_t_1.RData")
										PCL_raw_class[PCL_t_1] <- 1 
										rm (PCL_t_1);
									}	
									PCL_t_2 <- dget ("PCL_t_2.RData")		
										PCL_raw_class[PCL_t_2] <- 2
									rm (PCL_t_2);
	print ("+++++++++++++++++++++++++++++++++++PART 4 +++++++++++++++++++++++++++++++++++")
	print(Sys.time())

									PCL_t_3 <- dget ("PCL_t_3.RData")
										PCL_raw_class[PCL_t_3] <- 3
									rm (PCL_t_3);
									PCL_t_4 <- dget ("PCL_t_4.RData")
										PCL_raw_class[PCL_t_4] <- 4
									rm (PCL_t_4);
									PCL_raw_class[is.na(PCL_raw_class)] <- 0
								# Make a 0,1 mask		
									PCL <- PCL_raw_class
									rm (PCL_raw_class);
									PCL[PCL!=0] <- 1
									dput (PCL, "PCL.RData")
							## Cloud work
								# Clumping (group unique clouds)
									Band_1 <- raster (INPUT_BANDS [1])
									PCL_ras <- setValues(Band_1,PCL)
									writeRaster (PCL_ras, "PCL_ras.tif", overwrite = TRUE)
									rm (PCL); rm (Band_1);
									PCL_clump <- clump (PCL_ras)
									PCL_clump_vec <- PCL_clump[]
									PCL_clump.ncell <- ncell(PCL_clump)
									rm (PCL_clump)
								## Add 8 pixels around each cloud pixel in a clump
									PCL <- dget("PCL.RData")
									Cloud_pix_row <- which (PCL == 1) # Cloud pixels
									rm (PCL)
									dput (Cloud_pix_row, "Cloud_pix_row.RData")
							# Adjacency (outside the edge of cloud)
								PCL_ras <- raster ("PCL_ras.tif")
								PCL_adj_1 <- adjacency (PCL_ras, Cloud_pix_row, 1 : PCL_clump.ncell, directions = 8)
								rm (PCL_ras)
								rm (Cloud_pix_row);
								PCL_adj_v1 <- vector (length = PCL_clump.ncell)			
								PCL_adj_v1[PCL_adj_1[,2]] <- 2
								rm (PCL_adj_1)
								Cloud_pix_row <- dget ("Cloud_pix_row.RData")
								PCL_adj_v1[Cloud_pix_row] <- 1
								dput(PCL_adj_v1,"PCL_adj_v1.RData")
								rm (Cloud_pix_row)
							# Adjacency (at the edge of a cloud)
								Outer_edge_cloud_pix <- which (PCL_adj_v1 == 2)
								rm (PCL_adj_v1);
								PCL_ras <- raster ("PCL_ras.tif")
								PCL_adj_2 <- adjacency (PCL_ras, fromCells = Outer_edge_cloud_pix,toCells = 1:PCL_clump.ncell, directions = 8)
								rm (PCL_ras)
								PCL_adj_v2 <- vector(length = PCL_clump.ncell)	
								Cloud_pix_row <- dget ("Cloud_pix_row.RData")
								Inner_edge <- which (PCL_adj_2[,2] %in% Cloud_pix_row)
								rm (Cloud_pix_row)
								PCL_adj_v2[PCL_adj_2[,2]] <- 3
								rm (PCL_adj_2)
								PCL_adj_v2[PCL_adj_v2 != 3] <- 0
								PCL_adj_v2_2 <- rep (NA, length (PCL_adj_v2))
								PCL_adj_v1 <- dget ("PCL_adj_v1.RData")
								PCL_adj_v2_2[PCL_adj_v1==1&PCL_adj_v2==3] <- 1
								rm (PCL_adj_v2); rm (PCL_adj_v1)
								PCL_adj_v2_2[is.na(PCL_adj_v2_2)==TRUE] <- 0
								Cloud_edge_pix <- PCL_adj_v2_2
								rm (PCL_adj_v2_2)
							# Re-clump without edges
								PCL <- dget("PCL.RData")
								PCL_edge <- PCL
								rm (PCL)
								PCL_edge[which (Cloud_edge_pix == 1)] <- 0
								rm (Cloud_edge_pix)
								# Can't really clump if only 1 pixel
									Band_1 <- raster (INPUT_BANDS [1]) # Band 1		
									PCL_ras_edge <- setValues (Band_1,PCL_edge)
									if (length(PCL_edge[PCL_edge==1]) > 1){
										PCL_clump_edge <- clump (PCL_ras_edge)
										PCL_clump_vec_edge <- PCL_clump_edge[]
											} else {
										PCL_clump_edge <- PCL_ras_edge
										PCL_clump_vec_edge <- PCL_edge
									}
									rm (PCL_edge); rm (Band_1)
									dput(PCL_clump_vec_edge,"PCL_clump_vec_edge.RData")	
						## EQ 22/23 Cloud base 
							BT <- dget( "Band_6_vec_na.RData") # Band6 and BT are the same thing, just change of nomenclature
							BT_split <- split (BT, PCL_clump_vec_edge) # Split BT values
							rm (BT)
							# Length of clouds
								PCL_clump_vec_edge.table <- table (PCL_clump_vec_edge)
								PCL_clump_vec_edge_table_df <- data.frame (VAL = names (PCL_clump_vec_edge.table), COU = as.numeric (PCL_clump_vec_edge.table))
								rm (PCL_clump_vec_edge.table); rm (PCL_clump_vec_edge);
								Tcloud_base <- rep (NA, nrow (PCL_clump_vec_edge_table_df)) # Blank vector
						# Cloud base BT
							for(a in 1 : length(Tcloud_base)){
								#print (paste("Cloud base",a,"of",length(Tcloud_base)))
								BT_split_iter <- BT_split[[a]]
								# Cloud_id <- PCL_clump_vec_edge_table_df[a,1] # Not used anywhere else. Don't know why here. (JH)
								R <- sqrt (PCL_clump_vec_edge_table_df[a,2] / (2 * pi))	
								if (R >= 8){ #C- Between 0 and 1 in R -C#
									quant.val <- (((R - 8)^2)) / R^2
									Tcloud_base[a] <- as.numeric (quantile (BT_split_iter,quant.val, na.rm = TRUE))
								}
								if (R < 8){
									Tcloud_base[a] <- min (BT_split_iter)
									Tcloud_base[a] <- ifelse (Tcloud_base[a] < 0, NA, Tcloud_base[a]) #c- read through more -c#
								}
							}
							dput(Tcloud_base,"Tcloud_base.RData")
							rm (PCL_clump_vec_edge_table_df); rm (BT_split_iter); rm (R);
	print ("+++++++++++++++++++++++++++++++++++PART 5 +++++++++++++++++++++++++++++++++++")
	print(Sys.time())

						## EQ 23: Tcloud_object > base == base
							BT_split_mod <- BT_split
							for(a in 1:length (Tcloud_base)){
								#print (paste("Cloud base EQ23",a,"of",length(Tcloud_base)))
								in_1 <- BT_split[[a]]
								in_2 <- Tcloud_base[a]
								in_1[in_1 > in_2] <- in_2
								BT_split_mod[[a]] <- in_1
								rm (in_1); rm (in_2);
							}
							rm (BT_split)
							PCL_clump_vec_edge <- dget ("PCL_clump_vec_edge.RData")
							BT_mod <- unsplit (BT_split_mod, PCL_clump_vec_edge) # Modified BT
							rm (BT_split_mod); rm (PCL_clump_vec_edge);
							dput(BT_mod,"BT_mod.RData")
						## Eq 20: potential snow layer
							modified_NDSI <- dget ("modified_NDSI.RData")
							BT <- dget( "Band_6_vec_na.RData") # Band6 and BT are the same thing, just change of nomenclature
							Band_4_vec_na <- dget( "Band_4_vec_na.RData")
							Band_2_vec_na <- dget( "Band_2_vec_na.RData")
							PSL <- ifelse (modified_NDSI > 0.15 & BT < 380 & Band_4_vec_na > 1100 & Band_2_vec_na > 1000, TRUE, FALSE)
							rm (modified_NDSI); rm (BT); rm (Band_4_vec_na); rm (Band_2_vec_na); 
							dput(PSL,"PSL.Rdata")
							rm (PSL)
						## Object-based cloud and cloud shadow match
							PCL <- dget("PCL.RData")
							revised_ptm <- 100 * (length (which (PCL == 1)) / length (PCL))
							rm (PCL);
							if (revised_ptm >= 90){ # Too many clouds.  Set entire layer to cloud.
								print("No shadow match because too much cloud")
								Band_1_vec_na <- dget( "Band_1_vec_na.RData")
								Band_1 <- raster (INPUT_BANDS [1])
								Fmask <- rep (1, length (Band_1_vec_na)) # Set vector values
								Fmask_ras <- setValues (Band_1, Fmask)
								rm (Band_1); rm (Band_1_vec_na);
								writeRaster (Fmask_ras, OUTPUT, overwrite = TRUE)
							}
							if (revised_ptm < 90){ # Not too many clouds
					# RSAGA
					# Remove spikes and dips
					# Band 4
						Band_1 <- raster (INPUT_BANDS [1])
						Band_4_vec_na <- dget("Band_4_vec_na.RData")
						Band_4_vec_ref <- setValues (Band_1, as.numeric (Band_4_vec_na))
						writeRaster(Band_4_vec_ref,paste(getwd(), "BAND_4_TEMP_filled.sgrd",sep="/"), overwrite = TRUE)
						rm (Band_1); rm (Band_4_vec_na); rm (Band_4_vec_ref);
						rsaga.fill.sinks(paste(getwd(), "BAND_4_TEMP_filled.sgrd",sep="/"), paste(getwd(), "BAND_4_TEMP_filled.sgrd",sep="/"), method = "planchon.darboux.2001")
					# Band 5
						Band_1 <- raster (INPUT_BANDS [1])
						Band_5_vec_na <- dget("Band_5_vec_na.RData")
						Band_5_vec_ref <- setValues (Band_1, as.numeric (Band_5_vec_na))
						writeRaster(Band_5_vec_ref,  paste(getwd(), "BAND_5_TEMP_filled.sgrd",sep="/"), overwrite = TRUE)
						rm (Band_1); rm (Band_5_vec_na); rm (Band_5_vec_ref);
						rsaga.fill.sinks(paste(getwd(), "BAND_5_TEMP_filled.sgrd",sep="/"), paste(getwd(), "BAND_5_TEMP_filled.sgrd",sep="/"), method = "planchon.darboux.2001")	
					# Re-format data	
						rsaga.sgrd.to.esri(in.sgrd=c("BAND_4_TEMP_filled.sgrd", "BAND_5_TEMP_filled.sgrd"), out.grids=c("BAND_4_TEMP_filled.tif", "BAND_5_TEMP_filled.tif"), out.path=getwd(), prec=1)			
					## Eq 19
						Band_4_filled_vec_na <- raster("BAND_4_TEMP_filled.tif")[]
						Band_5_filled_vec_na <- raster("BAND_5_TEMP_filled.tif")[]
						Band_4_vec_na <- dget("Band_4_vec_na.RData")
						Band_5_vec_na <- dget("Band_5_vec_na.RData")
						
						shadow_prob <- apply (cbind (Band_4_filled_vec_na - Band_4_vec_na,Band_5_filled_vec_na - Band_5_vec_na), 1, function(x) min (x, na.rm = TRUE))
						rm (Band_4_filled_vec_na); rm (Band_5_filled_vec_na);
	#START###################################################
		# COMMENT
		# JH (28/9/12): Can't find PCSL used anywhere is fmask. Has to go somehwere?	
						PCSL <- shadow_prob > 200
	#END#####################################################					
						rm (shadow_prob);
						rm (PCSL);
						## Constants
							Tcloud_base <- dget ("Tcloud_base.RData")
							t_obj <- Tcloud_base
							dput(t_obj,"t_obj.RData")
							Tsimilar <- 0.3
							Tbuffer <- 0.98 # Treshold for matching buffering					
							num_cldoj <- 25 # Minimum matched cloud object (pixels)
							num_pix <- 8 # Number of inward pixels (240 M) for cloud base temperature
							rate_elapse <- 6.5 # Lapse rates
							rate_dlapse <- 9.8 # Lapse rates
							Max_cl_height <- 12000 	# Max cloud base height (m)
							Min_cl_height <- 200 # Min cloud base height (m)				
						## Cloud heights
							t_tempL <- dget ("t_tempL.RData")
							Min_cl_height_1 <- unlist (lapply (t_obj, function (x) max (Min_cl_height, 10 * (t_tempL - 400 - x) / rate_dlapse, na.rm = TRUE)))
							rm (t_tempL)
							Min_cl_height_1 <- ifelse (Min_cl_height_1 < 0, NA, Min_cl_height_1)
							t_tempH <- dget ("t_tempH.RData")
							Max_cl_height_1 <- unlist (lapply (t_obj, function (x) min (Max_cl_height, 10 * (t_tempH - 400 - x), na.rm = TRUE)))
							rm (t_tempH)
							Max_cl_height_1 <- ifelse (Max_cl_height_1 < 0, NA, Max_cl_height_1)
						## Object-based cloud and cloud shadow match
							BT_mod <- dget("BT_mod.RData")
							PCL_clump_vec_edge <- dget ("PCL_clump_vec_edge.RData")
							BT_mod_split <- split (BT_mod, PCL_clump_vec_edge) #C- Iteration for each cloud -C#
							rm (BT_mod)
	print ("+++++++++++++++++++++++++++++++++++PART 6 (may take a while) +++++++++++++++++++++++++++++++++++")
	print(Sys.time())


							PCL_clump_vec_edge.mem <- data.frame (ID = 1:length (PCL_clump_vec_edge), PCL_clump_vec_edge = PCL_clump_vec_edge) #C- made due to memory problems -C#
							rm (PCL_clump_vec_edge)
							shadow_cells_list <- list()	
							h <- list()
							YX_mod_out_list <- list()
							## Each cloud
							Band_1 <-  raster (INPUT_BANDS [1])
							sub_size <- xres (Band_1) # Resolution (assumes xres==yres)
							rm (Band_1)
							
				
		
			
	#START###################################################
		# COMMENT
		# JH (28/9/12): Parallel computing may speed up this loop. It goes through each cloud, which means that iterations (loop of a) are independent.
		#	Error in YX_mod[[1]][, ] : incorrect number of dimensions
			
			Tcloud_base <- dget ("Tcloud_base.RData")
			t_obj <- dget("t_obj.RData")
			for(a in 1:length (Tcloud_base)){
			print (paste(a,"of",length (Tcloud_base)))
				#print(paste(a,"of",length (Tcloud_base)))
				if (is.na (Min_cl_height_1[a]) == FALSE & is.na (Max_cl_height_1[a]) == FALSE){
					#shadow_cells_list_out_b <- test.func.c(a)
	# #START###################################################
	# # COMMENT
	# # JH (28/9/12): Problem--minimum cloud height is bigger than max height)	
					Min_cl_height_1[a] <- ifelse(Min_cl_height_1[a]>Max_cl_height_1[a],200,Min_cl_height_1) # Temporary work-around
					 base_h <- seq (Min_cl_height_1[a], Max_cl_height_1[a], 10)
	# #END#####################################################
			temp_obj <- BT_mod_split[[a]]
			t_obj_a <- t_obj[a]
			temp_obj_length <- length (temp_obj)
			h.raw <- 10 * (t_obj_a - temp_obj)	
			i_xy.raw <-  (sub_size * tan (L7.SUN.ELEVATION.TAN.RADIANS))
			X_move.raw 	<- cos (L7.SUN.AZIMUTH.TAN.RADIANS)
			Y_move.raw <-  sin (L7.SUN.AZIMUTH.TAN.RADIANS)
			#b.list.1 <- vector("list", length(base_h))
			b.list.1 <- sapply(base_h, function (x) {
					#print(paste(which(base_h == x), "of",length(base_h)))
					# Each height in the cloud
						h <- h.raw / rate_elapse + x
					return (h)
			})
			b.vec.1 <- b.list.1
			i_xy <- b.vec.1 / i_xy.raw
			# Movement (pixels)
				X_move <- (i_xy * X_move.raw) / 30
				Y_move <- (i_xy * Y_move.raw) / 30
			# Indentify cloud cells
				PCL_clump_vec_edge.mem.a <- PCL_clump_vec_edge.mem[,1]
				PCL_clump_vec_edge.mem.b <- PCL_clump_vec_edge.mem[,2]
				PCL_clump_vec_edge.mem.c <- which (PCL_clump_vec_edge.mem.b == a)
				PCL_clump_vec_edge.mem.d <- PCL_clump_vec_edge.mem.a[PCL_clump_vec_edge.mem.c]
				rm (PCL_clump_vec_edge.mem.a); rm (PCL_clump_vec_edge.mem.b);rm (PCL_clump_vec_edge.mem.c);
				cloud.cell.id <- as.numeric (na.omit (PCL_clump_vec_edge.mem.d))
				YX_ori <- rowColFromCell (PCL_clump_edge, cloud.cell.id)	
			# Identify shadow rowCol positions
				YX_mod <- vector("list", 2)
				length (YX_ori[, 1] - Y_move)
			if (L7.SUN.AZIMUTH < 180){
				YX_mod[[1]] <- YX_ori[,1] - Y_move
				YX_mod[[2]] <- YX_ori[,2] - X_move
				} else {
				YX_mod[[1]] <- YX_ori[,1] + Y_move
				YX_mod[[2]] <- YX_ori[,2] + X_move
				rm(YX_ori)
			}
			# Values that are off the map
				#YX_mod[YX_mod[,1] <= 0,] <- NA
				if (temp_obj_length == 1){
						YX_mod[[1]][YX_mod[[1]] <= 0] <- NA
						YX_mod[[2]][YX_mod[[2]] <= 0] <- NA
					} else {
						YX_mod[[1]][YX_mod[[1]][,] <= 0] <- NA
						YX_mod[[2]][YX_mod[[2]][,] <= 0] <- NA
					}
				#YX_mod[[1]] <- apply(YX_mod[[1]],2,function (x) ifelse (x <0, NA, x))
				#YX_mod[[2]] <- apply(YX_mod[[2]],2,function (x) ifelse (x <0, NA, x))
			# Merge YX_mod
				if (temp_obj_length == 1){
					o1 <- c(YX_mod[[1]],YX_mod[[2]])
					} else {
					o1 <- c(YX_mod[[1]][,1],YX_mod[[2]][,1])
				}
					o2 <- matrix (o1,ncol=2)
					YX_mod_merge_1 <- na.omit(o2)
					attributes(YX_mod_merge_1)$na.action <- NULL
					rm (o1); rm (o2)
				if (temp_obj_length != 1){	
					YX_mod_merge_2 <- vector("list", (ncol(YX_mod[[1]])-1))
						for(i in 1:(ncol(YX_mod[[1]])-1)){
							o1 <- c(YX_mod[[1]][,(i+1)],YX_mod[[2]][,(i+1)])
							o2 <- matrix (o1,ncol=2)
							rm (o1)
							YX_mod_merge_2[[i]] <- na.omit(o2)
							attributes(YX_mod_merge_2[[i]])$na.action <- NULL
						}
				}
			# Shadow cells	
				shadow_cells_list_out_b <- shad_cell_func (PCL_clump_edge, as.data.frame (YX_mod_merge_1, stringsAsFactors = FALSE))
				#count <- 0
				if (temp_obj_length != 1){	
					shadow_cells_list_out_b <-sapply(YX_mod_merge_2, function (x){
						#count <- ifelse(count!=0,count+1,1)
						#print (count)
						o1 <-  as.numeric (unique (na.omit (c (shadow_cells_list_out_b, shad_cell_func (PCL_clump_edge, as.data.frame (x, stringsAsFactors = FALSE))))))
						return (o1)
						})
				} else {
					shadow_cells_list_out_b <- as.numeric (unique (na.omit (shadow_cells_list_out_b)))
				}
				} else {
					shadow_cells_list_out_b <- NA
				}
			}
		# JH (29/9/12): Mistmatch for some vector lenghts	
	#END#####################################################
			dput(shadow_cells_list_out_b,"shadow_cells_list_out_b.RData")
	print ("+++++++++++++++++++++++++++++++++++PART 7 +++++++++++++++++++++++++++++++++++")
	print(Sys.time())

			## Cloud raster
				Unique.cloud.pix <- unique(shadow_cells_list) # Unique pixels only
				if (length(Unique.cloud.pix)!=0){
					cloud_shadow <- rep (NA, length (PCL_clump_edge[])) # Make vector
					cloud_shadow[unlist (Unique.cloud.pix)] <- 1 # Shadow	
					Band_1 <-  raster (INPUT_BANDS [1])						
					cloud_shadow_ras_1 <- setValues (Band_1, cloud_shadow) # Make raster
					rm (cloud_shadow)
					## Dialte cloud and shadow pixels 3 times. Zhu and Woodcock, 2012 call for buffers by 3 pixels in 8-connected directions -C#
					# 1 pixel out
						shadow_cal_1 <- adjacency (cloud_shadow_ras_1, fromCells = Unique.cloud.pix, toCells = c(1:ncell (cloud_shadow_ras_1)), directions = 8)
						rm (cloud_shadow_ras_1)
						shadow_cal_1_vec <- rep (NA, length (PCL_clump_edge[]))		
						shadow_cal_1_vec[shadow_cal_1[,2]] <- 1
						shadow_cal_1_ras <- setValues (Band_1, shadow_cal_1_vec)				
					# 2 pixels out
						shadow_cal_2_vec_vals <- which (shadow_cal_1_vec == 1)
						rm (shadow_cal_1_vec)
						shadow_cal_2 <- adjacency (shadow_cal_1_ras, fromCells = shadow_cal_2_vec_vals, toCells = c(1:ncell (shadow_cal_1)), directions = 8)
						rm (shadow_cal_1_ras)
						shadow_cal_2_vec <- rep (NA, length (PCL_clump_edge[]))		
						shadow_cal_2_vec[shadow_cal_2[,2]] <- 1
						rm (shadow_cal_2)
						shadow_cal_2_ras <- setValues (Band_1, shadow_cal_2_vec)
					# 3 pixels out
						shadow_cal_3_vec_vals <- which (shadow_cal_2_vec == 1)
						rm (shadow_cal_2_vec)
						shadow_cal_3 <- adjacency (shadow_cal_2_ras, fromCells = shadow_cal_3_vec_vals, toCells = c(1:ncell (shadow_cal_1)), directions = 8)
						rm (shadow_cal_2_ras)
						rm(shadow_cal_1)
						shadow_cal_3_vec <- rep (NA, length (PCL_clump_edge[]))		
						shadow_cal_3_vec[shadow_cal_3[,2]] <- 1
					} else {
						shadow_cal_3_vec <- rep (0, length (PCL_clump_edge[]))
					}
					
				## Classify output #C- Priority: 1) cloud pixels, 2) shadow pixels, 3) snow pixels
					Fmask <- rep(0,length(PCL_clump_edge[]))						
					PSL <- dget("PSL.RData")
					Fmask[PSL==TRUE] <- 3 # Snow
					rm (PSL)
					Fmask[shadow_cal_3_vec[]==1] <- 2 # Shadows
					rm (shadow_cal_3_vec)
					PCL <- dget("PCL.RData")
					Fmask[PCL==1] <- 1 # Cloud
					rm (PCL)
					Fmask[id_missing] <- NA # NA post 2003
					Band_1 <-  raster (INPUT_BANDS [1])	
					Fmask_ras <- setValues(Band_1,Fmask) # Raster
					writeRaster(Fmask_ras, OUTPUT, overwrite = TRUE)
				}
			}
		 }
	}
}

